Load packages
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Check working directory
getwd()
## [1] "/Users/jillashey/Desktop/PutnamLab/Repositories/Astrangia_repo/Astrangia_repo"
Load data
# Load data
daily <- read.csv("data/DailyMeasurements.csv", header = T)
daily <- daily[-13,] # remove tank 14 because it cracked
daily$Timepoint <- as.Date(daily$Timepoint, "%m/%d/%Y")
daily$Tank <- as.factor(daily$Tank)
daily <- subset(daily, Treatment=="Heat" | Treatment=="Ambient")
Discrete pH calculations from Tris calibrations.
path <-("data/pH_Tris") #set path to calibration files
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get the csv files
pH.cals <- data.frame(matrix(NA, nrow=length(file.names), ncol=3, dimnames=list(file.names,c("Date", "Intercept", "Slope")))) #generate a 3 column dataframe with specific column names
for(i in 1:length(file.names)) { # for every file in list start at the first and run this following function
Calib.Data <-read.table(file.path(path,file.names[i]), header=TRUE, sep=",", na.string="NA", as.is=TRUE) #reads in the data files
file.names[i]
model <-lm(mVTris ~ TTris, data=Calib.Data) #runs a linear regression of mV as a function of temperature
coe <- coef(model) #extracts the coeffecients
summary(model)$r.squared #extracts the r squared
plot(Calib.Data$mVTris, Calib.Data$TTris) #plots the regression data
pH.cals[i,2:3] <- coe #inserts coefficients in the dataframe
pH.cals[i,1] <- substr(file.names[i],1,8) #stores the file name in the Date column
}
colnames(pH.cals) <- c("Calib.Date", "Intercept", "Slope") #rename columns
pH.cals #view data
## Calib.Date Intercept Slope
## 20210122.csv 20210122 -86.73848 1.0218902
## 20210126.csv 20210126 -86.23930 0.9322418
## 20210203.csv 20210203 -88.14064 1.0017943
## 20210205.csv 20210205 -88.23311 0.8416322
## 20210217.csv 20210217 -87.53924 1.1237563
## 20210221.csv 20210221 -85.87848 0.8968835
## 20210223.csv 20210223 -86.36221 0.9370372
## 20210225.csv 20210225 -85.86760 0.9234509
## 20210227.csv 20210227 -87.60427 1.0150785
## 20210308.csv 20210308 -89.42691 0.9977500
## 20210317.csv 20210317 -86.21007 0.9799588
## 20210325.csv 20210325 -87.82550 1.0238389
## 20210409.csv 20210409 -88.31599 1.1153461
## 20210411.csv 20210411 -53.62103 -0.1385311
## 20210420.csv 20210420 -68.46091 0.2421390
## 20210524.csv 20210524 -96.95526 1.1877897
## 20210527.csv 20210527 -98.21886 1.2378231
## 20210531.csv 20210531 -88.25076 0.8102227
## 20210603.csv 20210603 -92.83213 1.0612981
## 20210609.csv 20210609 -97.02534 1.2354378
## 20210611.csv 20210611 -97.84557 1.2647915
## 20210615.csv 20210615 -96.14596 1.1781620
## 20210618.csv 20210618 -98.24662 1.2787030
## 20210622.csv 20210622 -96.86842 1.1817890
## 20210625.csv 20210625 -100.50568 1.3608393
## 20210628.csv 20210628 -92.85454 1.0585358
## 20210702.csv 20210702 -95.97090 1.1559734
## 20210706.csv 20210706 -101.52648 1.3059321
## 20210708.csv 20210708 -98.39602 1.1416130
## 20210713.csv 20210713 -102.79400 1.3403794
## 20210720.csv 20210720 -100.01209 1.0809837
#constants for use in pH calculation
R <- 8.31447215 #gas constant in J mol-1 K-1
F <- 96485.339924 #Faraday constant in coulombs mol-1
#merge with Seawater chemistry file
SW.chem <- merge(daily, pH.cals, by="Calib.Date")
Calculate total pH.
mvTris <- SW.chem$Temperature*SW.chem$Slope+SW.chem$Intercept #calculate the mV of the tris standard using the temperature mv relationships in the measured standard curves
STris<-35 #salinity of the Tris
phTris<- (11911.08-18.2499*STris-0.039336*STris^2)*(1/(SW.chem$Temperature+273.15))-366.27059+ 0.53993607*STris+0.00016329*STris^2+(64.52243-0.084041*STris)*log(SW.chem$Temperature+273.15)-0.11149858*(SW.chem$Temperature+273.15) #calculate the pH of the tris (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
SW.chem$pH.Total<-phTris+(mvTris/1000-SW.chem$pH.MV/1000)/(R*(SW.chem$Temperature+273.15)*log(10)/F) #calculate the pH on the total scale (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
Plot daily measurements
# By Treatment
temp.trt = ggplot(SW.chem, aes(x=Treatment, y=Temperature.C)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Temperature°C") +
theme(axis.text.x = element_text(angle = 90))
temp.trt
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
sal.trt = ggplot(SW.chem, aes(x=Treatment, y=Salinity.psu)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Salinity.psu") +
theme(axis.text.x = element_text(angle = 90))
sal.trt
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
pH.trt = ggplot(SW.chem, aes(x=Treatment, y=pH.Total)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("pH Total Scale") +
theme(axis.text.x = element_text(angle = 90))
pH.trt
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
light.trt = ggplot(SW.chem, aes(x=Treatment, y=Light)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Light µmol photons m^2 s^-1") +
theme(axis.text.x = element_text(angle = 90))
light.trt
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
flow.trt = ggplot(SW.chem, aes(x=Treatment, y=Flow.mL.sec)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Flow mL/sec") +
theme(axis.text.x = element_text(angle = 90))
flow.trt
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).
plot.trt <- grid.arrange(temp.trt, sal.trt, pH.trt, light.trt, flow.trt, ncol=3, nrow = 2, clip="off")
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).
ggsave("output/Daily_Measurements.By.Treatment.pdf", plot.trt, width = 21, height = 21, units = c("in"))
# By Tank
temp.tank = ggplot(SW.chem, aes(x=Tank, y=Temperature.C)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Temperature°C") +
theme(axis.text.x = element_text(angle = 90))
temp.tank
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
sal.tank = ggplot(SW.chem, aes(x=Tank, y=Salinity.psu)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Salinity.psu") +
theme(axis.text.x = element_text(angle = 90))
sal.tank
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
pH.tank = ggplot(SW.chem, aes(x=Tank, y=pH.Total)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("pH Total Scale") +
theme(axis.text.x = element_text(angle = 90))
pH.tank
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
light.tank = ggplot(SW.chem, aes(x=Tank, y=Light)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Light µmol photons m^2 s^-1") +
theme(axis.text.x = element_text(angle = 90))
light.tank
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
flow.tank = ggplot(SW.chem, aes(x=Tank, y=Flow.mL.sec)) +
geom_boxplot(aes(color = Treatment), size = 1) +
ylab("Flow mL/sec") +
theme(axis.text.x = element_text(angle = 90))
flow.tank
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).
plot.tank <- grid.arrange(temp.tank, sal.tank, pH.tank, light.tank, flow.tank, ncol=3, nrow = 2, clip="off")
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).
ggsave("output/Daily_Measurements.By.Tank.pdf", plot.tank, width = 21, height = 21, units = c("in"))
test for differences between tanks and treatments
# By Treatment
mod.temp <- aov(Temperature.C ~ Treatment, data = SW.chem)
mod.temp
## Call:
## aov(formula = Temperature.C ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 6140.63 53132.89
## Deg. of Freedom 1 1827
##
## Residual standard error: 5.392777
## Estimated effects may be unbalanced
## 11 observations deleted due to missingness
summary(mod.temp)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 6141 6141 211.1 <2e-16 ***
## Residuals 1827 53133 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness
hist(mod.temp$residuals)
mod.sal <- aov(Salinity.psu ~ Treatment, data = SW.chem)
mod.sal
## Call:
## aov(formula = Salinity.psu ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 17.6973 2115.0465
## Deg. of Freedom 1 1815
##
## Residual standard error: 1.079498
## Estimated effects may be unbalanced
## 23 observations deleted due to missingness
summary(mod.sal)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 17.7 17.697 15.19 0.000101 ***
## Residuals 1815 2115.0 1.165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23 observations deleted due to missingness
hist(mod.sal$residuals)
mod.pH <- aov(pH.Total ~ Treatment, data = SW.chem)
mod.pH
## Call:
## aov(formula = pH.Total ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 4.32348 91.46740
## Deg. of Freedom 1 1814
##
## Residual standard error: 0.2245508
## Estimated effects may be unbalanced
## 24 observations deleted due to missingness
summary(mod.pH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 4.32 4.323 85.74 <2e-16 ***
## Residuals 1814 91.47 0.050
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness
hist(mod.pH$residuals)
mod.light <- aov(Light ~ Treatment, data = SW.chem)
mod.light
## Call:
## aov(formula = Light ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 0.6 380831.6
## Deg. of Freedom 1 779
##
## Residual standard error: 22.11046
## Estimated effects may be unbalanced
## 1059 observations deleted due to missingness
summary(mod.light)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 1 0.6 0.001 0.972
## Residuals 779 380832 488.9
## 1059 observations deleted due to missingness
hist(mod.light$residuals)
mod.flow <- aov(Flow.mL.sec ~ Treatment, data = SW.chem)
mod.flow
## Call:
## aov(formula = Flow.mL.sec ~ Treatment, data = SW.chem)
##
## Terms:
## Treatment Residuals
## Sum of Squares 3167.9 679148.5
## Deg. of Freedom 1 78
##
## Residual standard error: 93.31148
## Estimated effects may be unbalanced
## 1760 observations deleted due to missingness
summary(mod.flow)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 3168 3168 0.364 0.548
## Residuals 78 679148 8707
## 1760 observations deleted due to missingness
hist(mod.flow$residuals)
# By Tank
mod.temp <- aov(Temperature.C ~ Tank, data = SW.chem)
mod.temp
## Call:
## aov(formula = Temperature.C ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 6154.91 53118.61
## Deg. of Freedom 15 1813
##
## Residual standard error: 5.412831
## Estimated effects may be unbalanced
## 11 observations deleted due to missingness
summary(mod.temp)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 6155 410.3 14.01 <2e-16 ***
## Residuals 1813 53119 29.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness
hist(mod.temp$residuals)
mod.sal <- aov(Salinity.psu ~ Tank, data = SW.chem)
mod.sal
## Call:
## aov(formula = Salinity.psu ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 18.8038 2113.9400
## Deg. of Freedom 15 1801
##
## Residual standard error: 1.083402
## Estimated effects may be unbalanced
## 23 observations deleted due to missingness
summary(mod.sal)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 18.8 1.254 1.068 0.382
## Residuals 1801 2113.9 1.174
## 23 observations deleted due to missingness
hist(mod.sal$residuals)
mod.pH <- aov(pH.Total ~ Tank, data = SW.chem)
mod.pH
## Call:
## aov(formula = pH.Total ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 4.57224 91.21864
## Deg. of Freedom 15 1800
##
## Residual standard error: 0.2251156
## Estimated effects may be unbalanced
## 24 observations deleted due to missingness
summary(mod.pH)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 4.57 0.30482 6.015 2.01e-12 ***
## Residuals 1800 91.22 0.05068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness
hist(mod.pH$residuals)
mod.light <- aov(Light ~ Tank, data = SW.chem)
mod.light
## Call:
## aov(formula = Light ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 17193.7 363638.5
## Deg. of Freedom 15 765
##
## Residual standard error: 21.80239
## Estimated effects may be unbalanced
## 1059 observations deleted due to missingness
summary(mod.light)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 17194 1146.2 2.411 0.00199 **
## Residuals 765 363638 475.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1059 observations deleted due to missingness
hist(mod.light$residuals)
mod.flow <- aov(Flow.mL.sec ~ Tank, data = SW.chem)
mod.flow
## Call:
## aov(formula = Flow.mL.sec ~ Tank, data = SW.chem)
##
## Terms:
## Tank Residuals
## Sum of Squares 51655.4 630661.0
## Deg. of Freedom 15 64
##
## Residual standard error: 99.26771
## Estimated effects may be unbalanced
## 1760 observations deleted due to missingness
summary(mod.flow)
## Df Sum Sq Mean Sq F value Pr(>F)
## Tank 15 51655 3444 0.349 0.987
## Residuals 64 630661 9854
## 1760 observations deleted due to missingness
hist(mod.flow$residuals)
Plot by date
# Plot and save
# By Treatment
temp.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
sal.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Salinity.psu, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
ph.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=pH.Total, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
light.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Light, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
flow.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Flow.mL.sec, group=Treatment)) +
geom_line(aes(color = Treatment), size = 1)
plot.trt.timepoint <- grid.arrange(temp.trt.date, sal.trt.date, ph.trt.date, light.trt.date, flow.trt.date, ncol=1, nrow = 5, clip="off")
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).
ggsave("output/Daily_Measurements.By.Date.Treatment.pdf", plot.trt.timepoint, width = 30, height = 21, units = c("in"))
# By Tank
temp.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Tank)) +
geom_line(aes(color = Tank), size = 1)
sal.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Salinity.psu, group=Tank)) +
geom_line(aes(color = Tank), size = 1)
ph.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=pH.Total, group=Tank)) +
geom_line(aes(color = Tank), size = 1)
light.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Light, group=Tank)) +
geom_line(aes(color = Tank), size = 1)
flow.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Flow.mL.sec, group=Tank)) +
geom_line(aes(color = Tank), size = 1)
plot.tank.timepoint <- grid.arrange(temp.tank.date, sal.tank.date, ph.tank.date, light.tank.date, flow.tank.date, ncol=1, nrow = 5, clip="off")
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).
ggsave("output/Daily_Measurements.By.Date.Tank.pdf", plot.tank.timepoint, width = 21, height = 25, units = c("in"))
plot.all <- grid.arrange(temp.trt.date, temp.tank.date, sal.trt.date, sal.tank.date, ph.trt.date,ph.tank.date, light.trt.date, light.tank.date, flow.trt.date, flow.tank.date, ncol=2, nrow=5)
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).
ggsave("output/Daily_Measurements.By.Date.pdf", plot.all, width = 30, height = 25, units = c("in"))